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Abstract 

A novel inversion technique is proposed to com- 
pute parametric maps showing the temperature, den- 
sity and chemical composition of cosmic hot gas from 
X-ray hyper-spectral images. The parameters are recov- 
ered by constructing a unique non-linear mapping derived 
by combining a physics-based modelling of the X-ray spec- 
trum with the selection of optimal bandpass filters. Prelim- 
inary results and analysis are presented. 

1. Introduction 

With the advent of space-borne X-ray astronomy in 
the 1960s, it became apparent that there are many cosmic 
sources emitting high energy electromagnetic radiation, in- 
dicating the presence of matter at very high temperatures. 
Point-like X-ray sources arise from hot gas falling onto col- 
lapsed objects: white dwarfs, neutron stars and black holes, 
but sources of diffuse X-ray emission were also discovered. 
Such extended sources generally indicate the presence of 
clouds of gas at temperatures above 10 6 K. Examples (on 
increasing size scales) include (a) supernova remnants, re- 
sulting from the explosion of massive dying stars, (b) hot 
gas within galaxies, often heated by supernova explosions, 
and (c) gas trapped in the deep gravitational potential wells 
of clusters of galaxies, typically 10 million light years in 
size, which are the largest stable structures in the Universe. 
We will deal with observations of (c) in this work. 

X-ray telescopes typically consist of a mirror system, fo- 
cusing X-rays onto a detector, which records the position 
and energy of each detected photon. The basic data product, 
when such a telescope is used to observe an extended cos- 
mic X-ray source, is a spectral map of the emission from the 
hot gas in the target. The spectral properties of this emission 
provide information about the temperature of the gas, and 
also about its chemical composition. In principle, therefore, 
information is available which would enable the tempera- 
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ture and composition of the gas to be mapped, enabling as- 
tronomers to investigate the origin of the gas, and the mech- 
anisms which have heated it. 

Unfortunately, unlike optical light for cosmic sources, X- 
ray photons are difficult to collect, since the earth's atmo- 
sphere absorbs and scatters them and thus X-ray observa- 
tories have to be space-borne. Furthermore, X-ray photons 
are scarce even from the brightest sources, and a long obser- 
vation of a comparatively bright source with a major X-ray 
observatory may detect only tens of thousands of photons. 
The issue of extracting reliable spatial-spectral information 
under these circumstances is not straightforward, and the 
methods commonly employed by astronomers have tended 
to be very CPU-intensive, and to produce less than satisfac- 
tory results, as will be described below. 

With the advent of large publicly available as- 
tronomical archives (the "Virtual Observatory", 
http://www.ivoa.net), it has become necessary to pro- 
duce algorithms for the fast and flexible extraction of maps 
of diffuse emission. We are developing a multi-spectral in- 
version approach which will allow such extraction of key 
gas parameters, complete with statistical errors, once an ini- 
tial inversion mapping has been calculated. A similar ap- 
proach could be employed in other applications, wherever 
a well-defined model is available, linking physical proper- 
ties to the resulting spectra. 

2. Physical mapping of hot gas in galaxy clus- 
ters 

The processes which generate X-ray emission in a hot 
plasma (a partially ionised gas) are complex and varied. Ex- 
tensive numerical codes (e.g. [1]) are available to calculate 
the intricate web of atomic excitations and de-excitations 
in the plasma, and compute the emergent spectrum un- 
der circumstances where the plasma is "optically thin" (re- 
absorption of photons in the hot plasma is negligible) and in 
thermodynamic equilibrium (i.e. has had sufficient time for 
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Figure 1 : (a) Adaptively smoothed X-ray image of the galaxy group centred on the galaxy NGC 5171, derived from data taken with the Eu- 
ropean Space Agency's XMM-Newton X-ray Observatory (http://xmm.vilspa.esa.es/); (b) Gas temperature map (white=hot to purple=cool) 
derived from the XMM-Newton observation of the galaxy group NGC5171 obtained by detailed spectral fitting. Contours show the sur- 
face brightness distribution from (a). 



the populations of the different ionisation states to achieve 
stability). Fortunately, these are good approximations in 
many astrophysical contexts, due to the very low density 
of the gas, and the long evolutionary timescales of the sys- 
tems concerned. 

2.1. Interpreting a spectral map of intergalactic 

gas 

In this paper we consider the task of mapping the phys- 
ical properties of the hot intergalactic gas within a small 
cluster (or "group") of galaxies. An X-ray image of this gas 
is shown in Fig. 1(a) for the group of galaxies surround- 
ing the nearby galaxy NGC 5171 [2]. The X-ray image 
shown in Fig. 1(a) has been processed in order to suppress 
the "Poisson noise", which arises from the limited num- 
ber of photons available, using an adaptive smoothing tech- 
nique. This involves convolving the raw image with a ker- 
nel whose width is small in regions of high brightness, but 
broad in regions of low brightness, where statistics are poor. 
The brightness of the emission is primarily an indication of 
the density of the hot gas, since the emissivity scales as the 
square of the gas density. 

In order to derive other physical parameters, it is neces- 
sary to extract spectra, and to compare these with one of the 
hot plasma models referred to above. These models are pa- 
rameterised, in simple cases, by a small number of param- 
eters - typically the gas temperature T, the abundance Z of 
heavy elements relative to the standard fractions seen within 
the solar system, and a normalisation constant from which 
the projected density of the gas can be derived. Absorption 
of some X-rays en route to us, by cool gas within our own 



Milky Way, acts to preferentially reduce the intensity of low 
energy X-rays, in a way which is well understood, and this 
introduces a further parameter into the spectral model (cor- 
responding to the projected density n# of cold gas along the 
line of sight to our target). Hence, in this simple case, the 
absorbed plasma model has four free parameters. For any 
choice of these parameters, the model can be used to pre- 
dict a spectrum which is incident on the telescope, this is 
then folded through the spectral response of the instrument, 
to produce a prediction for the detected spectrum, which 
is compared with that actually observed (after the removal 
of any background emission). The parameters of the spec- 
tral model are then optimised in an iterative process, and er- 
rors on these parameters derived from the way in which the 
fit statistic deteriorates as parameters are moved away from 
their best fit values. 

2.2. Current mapping method 

In order to map the physical properties of the gas over 
the image, it is normally necessary to extract spectra from 
a set of spatial regions [2] . However a minimum number of 
detected photons (approx. 1000) are required in order to de- 
rive meaningful constraints from the spectral fitting process. 
One therefore has to extract spectra from regions of varying 
size (larger in areas of low brightness), fit each individu- 
ally, and then assemble the results into an image in order 
to visualise the physical structure of the hot gas. This pro- 
cess is very time consuming, and produces a result which 
is very pixellated and hard to interpret (compare, for exam- 
ple, the blocky temperature map of Fig. 1(b) with the con- 
tinuous brightness distribution of Fig. 1(a)). 



3. A new parameter recovery method 



In a new method of parameter recovery, a set of corre- 
spondences is established between the physical parameters 
and multi-band image vectors. This is done in two steps. 
First, spectra are computed from parameters specifying ob- 
ject properties (see 2.1). Second, the equivalent image vec- 
tors are computed by convolving the spectra with suitably 
chosen bandpass filters. These two steps are computation- 
ally intensive, however they need to be carried out only once 
for a given imaging set-up. The set of correspondences es- 
tablished in this way forms the imaging model. The model 
is used to perform the inversion process, i.e. to infer the 
combination of parameters which lead to a particular im- 
age vector. This last step, which effectively performs image 
interpretation, is very fast and in some instances can be per- 
formed in real time. Its result is a set of parametric maps 
which show the magnitude of each parameter in a sepa- 
rate 2-dimensional image. In addition to improved speed, 
the ability to smooth the bandpass filtered images in a man- 
ner similar to that used in Fig. 1(a), will enable us to gen- 
erate continuously varying parameter maps representing far 
better the continuous distributions expected in the hot gas. 

3.1. Modelling of the spectra 

Using widely available software [1] synthetic spectra 
were generated for X-ray emitting optically thin plasmas 
with the range of physical parameters covering the range 
of interest, namely: the temperature T (0.6-1.5 keV), ele- 
mental abundance Z (0-0.9 times their corresponding solar 
values), and line-of-sight absorption due to a column of hy- 
drogen gas n H (0-1.8 x 10 25 m~ 2 ). Fig. 2 shows three ex- 
amples illustrating how the shape and detailed structure of 
the spectra in the range 0.2^- keV vary with various val- 
ues of these parameters. The fourth parameter, the overall 
normalisation, is adjusted to keep the total number of pho- 
tons detected to a constant value. The spectra were then con- 
volved with the telescope and detector response functions 
of the Advanced CCD Imaging spectrometer of the Chan- 
dra X-ray Observatory, the instrument used to obtain the ob- 
served data of HCG 62 used in the following section. 

3.2. Computational forward model 

Let an instance of hyper-spectral X-ray image data be 
represented by a discrete vector 

X = {X m },m = \,...,M, AeA (1) 

and the parameters depicting the physical properties be rep- 
resented by a vector 

p = peP (2) 
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Figure 2: Three examples of the spectra used to construct our fil- 
ters. They correspond to the following values for the three param- 
eters (T in keV, absorbing column n# in 10 25 m~ 2 , abundance Z 
as fraction of solar)= full line (0.6, 1.0, 0.3); dashed line (1.5, 0.4, 
0.6) and dotted line (1.2, 0.4, 0.3). 



The mapping / 

f-.P^A (3) 

models the physical processes governing the formation of 
spectra, for example the hot plasma model described in 2. 

In the first step of the construction of the imaging model, 
/ is used to generate a set of spectra corresponding to the 
entire range of the physically valid instances of parameters 
p, suitably discretised. Fig. 2 shows three examples of spec- 
tra generated using the model. 

In the existing parameter recovery method (see 2.2) the 
parameter vector p is recovered directly from the hyper- 
spectral data by iteratively minimising the distance between 
the measured and the computed spectrum. Given the com- 
plexity of the model / and the large dimensionality (M) of 
the spectra space A, this process is inevitably slow. Reduc- 
ing M is likely to speed-up the process, but at the cost of loss 
of accuracy. However, provided certain conditions are ful- 
filled, it is possible to reduce the dimensionality in such a 
way that relatively accurate information regarding the orig- 
inal parameter values can still be recovered [4]. 

Mapping g 

g : A — ► I IeR n (4) 

defines a projection from M-dimensional X-ray spectra 
space A into N-dimensional bandpass image space I such 
that N«M. The mapping is implemented as bandpass fil- 
tering, through a convolution of the spectral vector X € R M 
with N filter response functions. The response function for 



the n-th filter can be denned as 

M 

in = J^ C m^m,i = {i n },n = l,...,N (5) 

m=l 

The mapping function g can be defined by specifying the el- 
ements of matrix C" n such as to meet the required optimisa- 
tion criteria, as will be discussed in the next section. 
With the introduction of g, the mapping 

h = fog :P— »/ (6) 

describes the correspondence between specific param- 
eters characterising the astronomical object and the 
low-dimensional image vector computed by convolv- 
ing the large-dimensional hyper-spectral image through 
a set of optimally chosen filters C£. Whereas / is deter- 
mined by physics, and in this sense "fixed", g can be altered 
to meet user-defined optimisation criteria. 

3.3. Optimal spectral filter selection 

Drawing on the concepts introduced in 3.2 above, the 
task of estimating the parameters which characterise an ob- 
ject is to find an inverse function 

hT l :I — >P (7) 

such that hr l is unique, p = h~ x (i) exists for all peP 
and the error of mapping from the image data to the param- 
eter data is at its minimum. 

The uniqueness of mapping is the necessary condition 
and is tested by examining the behaviour of the determi- 
nant of the Jacobian matrix over all peP. One-to-one cor- 
respondence is ensured if and only if h is monotonic on the 
entire domain P of p. If h were continuous, this can be eas- 
ily determined by ensuring that all partial derivatives of h 
are strictly larger than zero over P. When, as in this case, h 
is a discrete vector valued function of a vector variable, the 
inverse function theorem [5] can be used to define an equiv- 
alent set of conditions: 

1. det{J)^Q for all peP, 

di n 

where J = — — , n=l,...,N; k=l,...,K, is a Jacobian ma- 

OPk 

trix of discrete partial derivatives 

2. the sign of det(J) is the same over the whole domain of 

P eP[4]. 

If these conditions are fulfilled, then, according to [5], 
there exists a neighbourhood around each pGP where there 
is one-to-one mapping between parameters p and image 
vectors i. Within such neighbourhood the inverse function 
h~ x can be expressed as 



where dp — p — po, di = i io and io = h(po). 

There could be many functions h meeting the unique- 
ness criterion. However, due to the errors inherent in both 
the model and the measurements, the parameter recovery 
errors will be different for different choices of h = fog. 
Since only the g component can be manipulated (see 3.2), 
the overall effect of h is determined by the choice of fil- 
ter functions, C^. The optimal mapping function is defined 
as the one which minimises the parameter recovery error. 

Using the appropriate statistical analysis [6] er- 
ror in the k-th parameter can be estimated to be: 



poiss\2\l c poiss_-^M rvi n 

Olfi 

P m is error corresponding to Poisson noise at energy level 

m, and — are available from the inverse of the Jacobian 
di n 

matrix, J~ l . The overall error is computed as the sum of % 
The outline of the algorithm used for finding opti- 
mal filters is given below. The optimisation procedure 
was implemented using a technique based on Evolution- 
ary Strategy [7], although other optimization algorithms 
could have been used with similar effect. Simple rect- 
angular filters were used, each filter defined by its mini- 
mum and maximum energy. 

Until the stopping criterion is met 

1 . define a new set of filters C n m 

2. for a given set of Q, and for each discrete p com- 
pute i = g(C£, f(p)) 

3. check that J is either strictly positive or strictly nega- 
tive for ALL p 

if true, compute the inverse Jacobian matrix 
if false, return to step 1 

4. compute the error of parameter recovery and use it to 
compute a stopping criterion. 

The use of the Jacobian matrix requires that the number 
of filters (N) is equal to the number of parameters (K). One 
of the parameters, the normalisation constant, is required to 
compensate for the uneven energy levels across the image 
space due to the nature of the imaging process and proper- 
ties of the imaging equipment. For this reason, instead of 
computing the image vector in step 2, an image quotient 

in 

vector, defined as q n = — , n=l,...N-l was computed and 
in 

the Jacobian matrix was defined as 

J=^-,n = l,...,N-\;k=l,...,K;K=N-\ (9) 



dp = J l di 



(8) 




Figure 3: A hyper-spectral X-ray image of the compact galaxy 
group HCG 62, in the constellation Virgo. The image is 4 arcmin 
on a side. Colours represent various levels of X-ray surface bright- 
ness: green depicts lower brightness, while purple indicates in- 
creasing X-ray intensity. 



3.4. Implementing mapping functions h and h 

The selection of optimal filters completes the modelling 
process. All the spectra generated with the model / and cor- 
responding to all the valid parameter vectors p (see 3.1) are 
now convolved with the filters (mapping g, see 3.3), the re- 
sulting image quotients q are then computed and stored in a 
discrete look-up table indexed by p. This table implements 
the (forward) mapping function h, such that q = /i(p) for all 
p £ P. As a set of parameter vectors p used to generate a 
specific embodiment of h is discrete, the computation of h 
for the parameters which are not in the generating set in- 
volves suitable interpolation (e.g. [8]). 

The inverse mapping function h : / — ► P is imple- 
mented as a simple look-up, i.e. given a quotient vector 
q = {j^} , n=l,...N-l„ find p^. which minimises the differ- 
ence ||q„— h(pt) ||. As in the forward case, interpolation will 
produce a more accurate result. 

3.5. Creating parametric maps 

Parametric maps are computed as follows: 

1 . Pre-process hyper-spectral images (compensate for read- 
out noise, remove instrument related artefacts etc.) 

2. Filter the images through to obtain N bandpass im- 
ages 

3. Apply spatial smoothing to the bandpass images 

4. Apply exposure correction 

5. Compute image quotients by dividing the images in the 
first N-l bands by the image in N th band 

6. For each image location recover the parameter vector 



from the image quotient, p = h~ x (q) 

7. Store each component of the parameter vector p,t, in a 
separate image (a parameteric map k) 
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Figure 4: A selection of spectra showing variation with metal- 
licity (Z=0 to 0.9 solar) with fixed T=0.6 keV and fixed njj = 
0.55 x 10 25 m~ 2 . Optimal filter set is shown overlaid. 



4. Results 

The inversion process outlined in the previous section 
was applied to a hyper-spectral X-ray image [9], shown in 
Fig. 3. The raw observation data was reduced using stan- 
dard techniques to correct for instrumental response, and 
the background detected from the sky and detector was 
subtracted using blank field observations. Spectra were ex- 
tracted in the same energy range as the synthetic spectra, 
from seven annuli, centred on the centroid of emission, of 
radii 20, 36, 52, 70, 1 10, 170 and 252 arcsec. 

4.1. Optimal filters 

In order to construct the optimal filters, 1000 synthetic 
spectra were generated for X-ray emitting optically thin 
plasmas as described in section 3.1. Parameter space was 
discretised into a 10x10x10 grid linearly spaced in each pa- 
rameter to cover the respective parameter ranges (T=0.6 to 
1.5 keV, Z=0.0 to 0.9 and n H =0 to 1.8 x I0 25 m - 2 ). Poisson 
errors, corresponding to a total count of 10 4 per spectrum, 
were assumed in order to define the wavelength dependent 
error associated with the spectral data. 

Following the process outlined in §3.3 and after approxi- 
mately 1 .3 x 10 7 Jacobian constructions, the filter set which 
both guarantees a non zero Jacobian and minimises the sum 
(over model parameter space) of parameter recovery errors 
was found and is shown in Fig. 4 overlaid on a selection of 
model spectra. The four filters in the set are marked Fl to 
F4 with F4 used to construct the denominator in quotient 
space. 



4.2. Parameter recovery 

In order to test the technique, parameter recovery was 
implemented for the spectra extracted from the seven an- 
nuli of the galaxy group HCG 62 (Fig. 3). The spectrum 
from each region was obtained by integrating the signal over 
that region. Using the optimal filter set derived in 4.1 the 
quotient vector for each region was produced and then, us- 
ing the inverse method (section 3.3), back projected to yield 
T, Z and n# and their associated errors. The results for the 
seven annuli are depicted in Figs. 5(a), 5(b) and 5(c) along 
with the results of the original algorithm for comparison. 
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Figure 5: Recovered parameter values and associated errors in the 
seven annular regions for (a) temperature (T); (b) metallicity (Z); 
and (c) intervening gas column (n#). 



5. Discussion and conclusions 

We have proposed a fast inversion technique which al- 
lows most of the essential physical information present in 
X-ray spectra of hot gas to be presented in the form of para- 
metric maps. This approach improves on the existing tech- 
nique in two aspects. First, the laborious task of iteratively 
adjusting the three parameters and a scaling factor to "best 
fit" a modelled spectrum to experimental data can be re- 
placed with a simple application of pre-calculated filters to 
the experimental spectrum and the inverse mapping of the 
subsequent quotient vector to parameter space. Secondly, 
since the images can be manipulated using standard image 
processing techniques before the inversion is performed, 
our approach permits us to produce stabilised and contin- 
uous maps of the physical properties which are important 
to astrophysicists. This provides a more useful product than 
current techniques, and involves less computational effort. 

At present the mapping procedure is not guaranteed to 
produce perfect results (Fig. 5) but the technique was orig- 



inally designed for applications in visible light where the 
use of contiguous normalized filters was necessary to en- 
able practical application in digital photography. Such limi- 
tations are not relevant to this work and we already envisage 
splitting filters to make best use of knowledge of the char- 
acteristic emission energies of important metallic heavy el- 
ements in the hot plasmas of interest. Further improvements 
in the inverse mapping can be achieved by choosing a vari- 
able separation in our parameter space distributions (when 
constructing the model) so as to more evenly populate quo- 
tient space. Both these straightforward modifications will 
only introduce delays into the one-off filter construction and 
mapping function definition and will have no impact on the 
time required to invert experimental spectra. 

The method is generic and has been successfully ap- 
plied to extract parameters of diagnostic importance from 
the colour images of the skin [4]. 
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